# coding=utf-8
# Copyright 2024 The Google Research Authors.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

# -*- coding: utf-8 -*-
"""so8_supergravity.ipynb

Automatically generated by Colaboratory.

Copyright 2020 Google LLC

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    https://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
"""

# Let us explicitly ask for TensorFlow2.
# This installs a lot of stuff - and will take a while.
# !pip3 install tensorflow==2.0.1

import base64
import collections
import dataclasses
import hashlib
import itertools
import math
import numpy
import pprint
import scipy.optimize
import sys

import tensorflow as tf


print('TF version is:', tf.__version__)
print('NumPy version is:', numpy.__version__)


@dataclasses.dataclass(frozen=True)
class Solution(object):
  potential: float
  stationarity: float
  pos: numpy.ndarray


def np_esum(spec, *arrays, optimize='greedy'):
  """Numpy einsum with default greedy optimization."""
  return numpy.einsum(spec, *arrays, optimize=optimize)


def get_onb_transform(k_ab):
  if not numpy.allclose(k_ab, k_ab.real) or not numpy.allclose(k_ab, k_ab.T):
    raise ValueError('Bad Gramian.')
  eigvals, eigvecsT = numpy.linalg.eigh(k_ab)
  if not all(v * eigvals[0] > 0 for v in eigvals):
    raise ValueError('Non-definite Gramian.')
  onb_transform = numpy.einsum('a,na->an', eigvals**(-.5), eigvecsT)
  g = np_esum('ab,Aa,Bb->AB', k_ab, onb_transform, onb_transform)
  assert numpy.allclose(
    g, numpy.eye(g.shape[0]) * ((-1, 1)[int(eigvals[0] > 0)])
    ), 'Bad ONB-transform.'
  return onb_transform, numpy.linalg.inv(onb_transform)


def numpy_signature(a, digits=3):
  """Produces a signature-fingerprint of a numpy array."""
  # Hack to ensure that -0.0 gets consistently shown as 0.0.
  minus_zero_hack = 1e-100+1e-100j
  return base64.b64encode(
      hashlib.sha256(
          str((a.shape,
               ','.join(repr(x)
               for x in numpy.round(a + minus_zero_hack, digits).flat))
          ).encode('utf-8')
      ).digest()).decode('utf-8').strip('\n=')


def tformat(array,
            name=None,
            elem_filter=lambda x: abs(x) > 1e-8,
            fmt='%s',
            max_rows=numpy.inf,
            cols=120):
    """Formats a numpy-array in human readable table form."""
    # Leading row will be replaced if caller asked for a name-row.
    dim_widths = [
        max(1, int(math.ceil(math.log(dim + 1e-100, 10))))
        for dim in array.shape]
    format_str = '%s: %s' % (' '.join('%%%dd' % w for w in dim_widths), fmt)
    rows = []
    for indices in itertools.product(*[range(dim) for dim in array.shape]):
        v = array[indices]
        if elem_filter(v):
            rows.append(format_str % (indices + (v, )))
    num_entries = len(rows)
    if num_entries > max_rows:
      rows = rows[:max_rows]
    if cols is not None:
      width = max(map(len, rows))
      num_cols = max(1, cols // (3 + width))
      num_xrows = int(math.ceil(len(rows) / num_cols))
      padded = [('%%-%ds' % width) % s
                for s in rows + [''] * (num_cols * num_xrows - len(rows))]
      table = numpy.array(padded, dtype=object).reshape(num_cols, num_xrows).T
      xrows = [' | '.join(row) for row in table]
    else:
      xrows = rows
    if name is not None:
      return '\n'.join(
          ['=== %s, shape=%r, %d%s / %d non-small entries ===' % (
              name, array.shape,
            num_entries,
            '' if num_entries == len(rows) else '(%d shown)' % num_entries,
            array.size)] +
          [r.strip() for r in xrows])
    return '\n'.join(xrows)


def tprint(array, sep=' ', end='\n', file=sys.stdout, **tformat_kwargs):
    """Prints a numpy array in human readable table form."""
    print(tformat(array, **tformat_kwargs), sep=sep, end=end, file=file)


### Lie Algebra definitions for Spin(8), SU(8), E7.

def permutation_sign(p):
  q = [x for x in p]  # Copy to list.
  parity = 1
  for n in range(len(p)):
    while n != q[n]:
      qn = q[n]
      q[n], q[qn] = q[qn], q[n]  # Swap to make q[qn] = qn.
      parity = -parity
  return parity


def asymm2(a, einsum_spec):
  """Antisymmetrizes."""
  return 0.5 * (a - numpy.einsum(einsum_spec, a))


class Spin8(object):
  """Container class for Spin(8) tensor invariants."""

  def __init__(self):
    r8 = range(8)
    self.gamma_vsc = gamma_vsc = self._get_gamma_vsc()
    #
    # The gamma^{ab}_{alpha beta} tensor that translates between antisymmetric
    # matrices over vectors [ij] and antisymmetric matrices over spinors [sS].
    self.gamma_vvss = asymm2(
        numpy.einsum('isc,jSc->ijsS', gamma_vsc, gamma_vsc), 'ijsS->jisS')
    # The gamma^{ab}_{alpha* beta*} tensor that translates between antisymmetric
    # matrices over vectors [ij] and antisymmetric matrices over cospinors [cC].
    self.gamma_vvcc = asymm2(
        numpy.einsum('isc,jsC->ijcC', gamma_vsc, gamma_vsc), 'ijcC->jicC')
    #
    # The gamma^{ijkl}_{alpha beta} tensor that translates between antisymmetric
    # 4-forms [ijkl] and symmetric traceless matrices over the spinors (sS).
    g_ijsS = numpy.einsum('isc,jSc->ijsS', self.gamma_vsc, self.gamma_vsc)
    g_ijcC = numpy.einsum('isc,jsC->ijcC', self.gamma_vsc, self.gamma_vsc)
    g_ijklsS = numpy.einsum('ijst,kltS->ijklsS', g_ijsS, g_ijsS)
    g_ijklcC = numpy.einsum('ijcd,kldC->ijklcC', g_ijcC, g_ijcC)
    gamma_vvvvss = numpy.zeros([8] * 6)
    gamma_vvvvcc = numpy.zeros([8] * 6)
    for perm in itertools.permutations(range(4)):
      perm_ijkl = ''.join('ijkl'[p] for p in perm)
      sign = permutation_sign(perm)
      gamma_vvvvss += sign * numpy.einsum(perm_ijkl + 'sS->ijklsS', g_ijklsS)
      gamma_vvvvcc += sign * numpy.einsum(perm_ijkl + 'cC->ijklcC', g_ijklcC)
    self.gamma_vvvvss = gamma_vvvvss / 24.0
    self.gamma_vvvvcc = gamma_vvvvcc / 24.0

  def _get_gamma_vsc(self):
    """Computes SO(8) gamma-matrices."""
    # Conventions match Green, Schwarz, Witten's, but with index-counting
    # starting at zero.
    entries = (
        "007+ 016- 025- 034+ 043- 052+ 061+ 070- "
        "101+ 110- 123- 132+ 145+ 154- 167- 176+ "
        "204+ 215- 226+ 237- 240- 251+ 262- 273+ "
        "302+ 313+ 320- 331- 346- 357- 364+ 375+ "
        "403+ 412- 421+ 430- 447+ 456- 465+ 474- "
        "505+ 514+ 527+ 536+ 541- 550- 563- 572- "
        "606+ 617+ 624- 635- 642+ 653+ 660- 671- "
        "700+ 711+ 722+ 733+ 744+ 755+ 766+ 777+")
    ret = numpy.zeros([8, 8, 8])
    for ijkc in entries.split():
     ijk = tuple(map(int, ijkc[:-1]))
     ret[ijk] = +1 if ijkc[-1] == '+' else -1
    return ret


class SU8(object):
  """Container class for su(8) tensor invariants."""

  def __init__(self):
    # Tensor that translates between adjoint indices 'a' and
    # (vector) x (vector) indices 'ij'
    ij_map = [(i, j) for i in range(8) for j in range(8) if i < j]
    #
    # We also need the mapping between 8 x 8 and 35 representations, using
    # common conventions for a basis of the 35-representation, and likewise
    # for 8 x 8 and 28.
    m_35_8_8 = numpy.zeros([35, 8, 8], dtype=numpy.complex128)
    m_28_8_8 = numpy.zeros([28, 8, 8], dtype=numpy.complex128)
    for n in range(7):
      m_35_8_8[n, n, n] = +1.0
      m_35_8_8[n, n + 1, n + 1] = -1.0
    for a, (m, n) in enumerate(ij_map):
      m_35_8_8[a + 7, m, n] = m_35_8_8[a + 7, n, m] = 1.0
      m_28_8_8[a, m, n] = 1.0
      m_28_8_8[a, n, m] = -1.0
    #
    # The su8 'Generator Matrices'.
    t_aij = numpy.zeros([63, 8, 8], dtype=numpy.complex128)
    t_aij[:35, :, :] = 1.0j * m_35_8_8
    for a, (i, j) in enumerate(ij_map):
      t_aij[a + 35, i, j] = -1.0
      t_aij[a + 35, j, i] = 1.0
    self.ij_map = ij_map
    self.m_35_8_8 = m_35_8_8
    self.m_28_8_8 = m_28_8_8
    self.t_aij = t_aij


class E7(object):
  """Container class for e7 tensor invariants."""

  def __init__(self, spin8, su8):
    self._spin8 = spin8
    self._su8 = su8
    ij_map = su8.ij_map
    t_a_ij_kl = numpy.zeros([133, 56, 56], dtype=numpy.complex128)
    t_a_ij_kl[:35, 28:, :28] = (1 / 8.0) * (
        np_esum('ijklsS,qsS,Iij,Kkl->qIK',
                spin8.gamma_vvvvss, su8.m_35_8_8, su8.m_28_8_8, su8.m_28_8_8))
    t_a_ij_kl[:35, :28, 28:] = t_a_ij_kl[:35, 28:, :28]
    t_a_ij_kl[35:70, 28:, :28] = (1.0j / 8.0) * (
        np_esum('ijklcC,qcC,Iij,Kkl->qIK',
                spin8.gamma_vvvvcc, su8.m_35_8_8, su8.m_28_8_8, su8.m_28_8_8))
    t_a_ij_kl[35:70, :28, 28:] = -t_a_ij_kl[35:70, 28:, :28]
    #
    # We need to find the action of the su(8) algebra on the
    # 28-representation.
    su8_28 = 2 * np_esum('aij,mn,Iim,Jjn->aIJ',
                          su8.t_aij,
                          numpy.eye(8, dtype=numpy.complex128),
                          su8.m_28_8_8, su8.m_28_8_8)
    t_a_ij_kl[70:, :28, :28] = su8_28
    t_a_ij_kl[70:, 28:, 28:] = su8_28.conjugate()
    self.t_a_ij_kl = t_a_ij_kl
    #
    self.k_ab = numpy.einsum('aMN,bNM->ab', t_a_ij_kl, t_a_ij_kl)
    self.v70_as_sc8x8 = numpy.einsum('sc,xab->sxcab',
                                      numpy.eye(2),
                                      su8.m_35_8_8).reshape(70, 2, 8, 8)
    # For e7, there actually is a better orthonormal basis:
    # the sd/asd 4-forms. The approach used here however readily generalizes
    # to all other groups.
    self.v70_onb_onbinv = get_onb_transform(self.k_ab[:70, :70])


def get_proj_35_8888(want_selfdual=True):
  """Computes the (35, 8, 8, 8, 8)-projector to the (anti)self-dual 4-forms."""
  # We first need some basis for the 35 self-dual 4-forms.
  # Our convention is that we lexicographically list those 8-choose-4
  # combinations that contain the index 0.
  sign_selfdual = 1 if want_selfdual else -1
  ret = numpy.zeros([35, 8, 8, 8, 8], dtype=numpy.float64)
  #
  def get_selfdual(ijkl):
    mnpq = tuple(n for n in range(8) if n not in ijkl)
    return (sign_selfdual * permutation_sign(ijkl + mnpq),
            ijkl, mnpq)
  selfduals = [get_selfdual(ijkl)
               for ijkl in itertools.combinations(range(8), 4)
               if 0 in ijkl]
  for num_sd, (sign_sd, ijkl, mnpq) in enumerate(selfduals):
    for abcd in itertools.permutations(range(4)):
      sign_abcd = permutation_sign(abcd)
      ret[num_sd,
          ijkl[abcd[0]],
          ijkl[abcd[1]],
          ijkl[abcd[2]],
          ijkl[abcd[3]]] = sign_abcd
      ret[num_sd,
          mnpq[abcd[0]],
          mnpq[abcd[1]],
          mnpq[abcd[2]],
          mnpq[abcd[3]]] = sign_abcd * sign_sd
  return ret / 24.0


spin8 = Spin8()
su8 = SU8()
e7 = E7(spin8, su8)


assert (numpy_signature(e7.t_a_ij_kl) ==
        'MMExYjC6Qr6gunZIYfRLLgM2PDtwUDYujBNzAIukAVY'), 'Bad E7(7) definitions.'

### SO(p, 8-p) gaugings

def get_so_pq_E(p=8):
  if p == 8 or p == 0:
    return numpy.eye(56, dtype=complex)
  q = 8 - p
  pq_ratio = p / q
  x88 = numpy.diag([-1.0] * p + [1.0 * pq_ratio] * q)
  t = 0.25j * numpy.pi / (1 + pq_ratio)
  k_ab = numpy.einsum('aij,bij->ab', su8.m_35_8_8, su8.m_35_8_8)
  v35 = numpy.einsum('mab,ab,mM->M', su8.m_35_8_8, x88, numpy.linalg.inv(k_ab))
  gen_E = numpy.einsum(
    'aMN,a->NM',
    e7.t_a_ij_kl,
    numpy.pad(v35, [(0, 133 - 35)], 'constant'))
  return scipy.linalg.expm(-t * gen_E)

### Supergravity.

@dataclasses.dataclass(frozen=True)
class SUGRATensors(object):
  v70: tf.Tensor
  vielbein: tf.Tensor
  tee_tensor: tf.Tensor
  a1: tf.Tensor
  a2: tf.Tensor
  potential: tf.Tensor


def get_tf_stationarity(fn_potential, **fn_potential_kwargs):
  """Returns a @tf.function that computes |grad potential|^2."""
  @tf.function
  def stationarity(pos):
    tape = tf.GradientTape()
    with tape:
      tape.watch(pos)
      potential = fn_potential(pos, **fn_potential_kwargs)
    grad_potential = tape.gradient(potential, pos)
    return tf.reduce_sum(grad_potential * grad_potential)
  return stationarity


@tf.function
def dwn_stationarity(t_a1, t_a2):
  """Computes the de Wit-Nicolai stationarity-condition tensor."""
  # See: https://arxiv.org/pdf/1302.6219.pdf, text after (3.2).
  t_x0 = (
      +4.0 * tf.einsum('mi,mjkl->ijkl', t_a1, t_a2)
      -3.0 * tf.einsum('mnij,nklm->ijkl', t_a2, t_a2))
  t_x0_real = tf.math.real(t_x0)
  t_x0_imag = tf.math.imag(t_x0)
  tc_sd = tf.constant(get_proj_35_8888(True))
  tc_asd = tf.constant(get_proj_35_8888(False))
  t_x_real_sd = tf.einsum('aijkl,ijkl->a', tc_sd, t_x0_real)
  t_x_imag_asd = tf.einsum('aijkl,ijkl->a', tc_asd, t_x0_imag)
  return (tf.einsum('a,a->', t_x_real_sd, t_x_real_sd) +
          tf.einsum('a,a->', t_x_imag_asd, t_x_imag_asd))


def tf_sugra_tensors(t_v70, compute_masses, t_lhs_vielbein, t_rhs_E):
  """Returns key tensors for D=4 supergravity."""
  tc_28_8_8 = tf.constant(su8.m_28_8_8)
  t_e7_generator_v70 = tf.einsum(
      'v,vIJ->JI',
      tf.complex(t_v70, tf.constant([0.0] * 70, dtype=tf.float64)),
      tf.constant(e7.t_a_ij_kl[:70, :, :], dtype=tf.complex128))
  t_complex_vielbein0 = tf.linalg.expm(t_e7_generator_v70) @ t_rhs_E
  if compute_masses:
    t_complex_vielbein = t_lhs_vielbein @ t_complex_vielbein0
  else:
    t_complex_vielbein = t_complex_vielbein0
  @tf.function
  def expand_ijkl(t_ab):
    return 0.5 * tf.einsum(
        'ijB,BIJ->ijIJ',
        tf.einsum('AB,Aij->ijB', t_ab, tc_28_8_8),
        tc_28_8_8)
  #
  t_u_ijIJ = expand_ijkl(t_complex_vielbein[:28, :28])
  t_u_klKL = expand_ijkl(t_complex_vielbein[28:, 28:])
  t_v_ijKL = expand_ijkl(t_complex_vielbein[:28, 28:])
  t_v_klIJ = expand_ijkl(t_complex_vielbein[28:, :28])
  #
  t_uv = t_u_klKL + t_v_klIJ
  t_uuvv = (tf.einsum('lmJK,kmKI->lkIJ', t_u_ijIJ, t_u_klKL) -
            tf.einsum('lmJK,kmKI->lkIJ', t_v_ijKL, t_v_klIJ))
  t_T = tf.einsum('ijIJ,lkIJ->lkij', t_uv, t_uuvv)
  t_A1 = (-4.0 / 21.0) * tf.linalg.trace(tf.einsum('mijn->ijmn', t_T))
  t_A2 = (-4.0 / (3 * 3)) * (
      # Antisymmetrize in last 3 indices, taking into account antisymmetry
      # in last two indices.
      t_T
      + tf.einsum('lijk->ljki', t_T)
      + tf.einsum('lijk->lkij', t_T))
  t_A1_real = tf.math.real(t_A1)
  t_A1_imag = tf.math.imag(t_A1)
  t_A2_real = tf.math.real(t_A2)
  t_A2_imag = tf.math.imag(t_A2)
  t_A1_potential = (-3.0 / 4) * (
      tf.einsum('ij,ij->', t_A1_real, t_A1_real) +
      tf.einsum('ij,ij->', t_A1_imag, t_A1_imag))
  t_A2_potential = (1.0 / 24) * (
      tf.einsum('ijkl,ijkl->', t_A2_real, t_A2_real) +
      tf.einsum('ijkl,ijkl->', t_A2_imag, t_A2_imag))
  t_potential = t_A1_potential + t_A2_potential
  #
  return t_v70, t_complex_vielbein, t_T, t_A1, t_A2, t_potential


def so8_sugra_tensors(t_v70, tc_rhs_E):
  t_v70, t_complex_vielbein, t_T, t_A1, t_A2, t_potential = (
     tf_sugra_tensors(t_v70, False, 0.0, tc_rhs_E))
  return SUGRATensors(
      v70=t_v70,
      vielbein=t_complex_vielbein,
      tee_tensor=t_T,
      a1=t_A1,
      a2=t_A2,
      potential=t_potential)


def so8_sugra_scalar_masses(v70, so_pq_p):
  # Note: In some situations, small deviations in the input give quite
  # noticeable deviations in the scalar mass-spectrum.
  # Getting reliable numbers here really requires satisfying
  # the stationarity-condition to high accuracy.
  tc_rhs_E = tf.constant(get_so_pq_E(so_pq_p), dtype=tf.complex128)
  tc_e7_onb = tf.constant(e7.v70_onb_onbinv[0], dtype=tf.complex128)
  tc_e7_taMN = tf.constant(e7.t_a_ij_kl[:70, :, :], dtype=tf.complex128)
  t_v70 = tf.constant(v70, dtype=tf.float64)
  #
  def tf_grad_potential_lhs_onb(t_d_v70_onb):
    tape = tf.GradientTape()
    with tape:
      tape.watch(t_d_v70_onb)
      t_d_gen_e7 = tf.einsum(
          'a,aMN->NM',
          tf.einsum('Aa,A->a',
                    tc_e7_onb,
                    tf.complex(t_d_v70_onb, tf.zeros_like(t_d_v70_onb))),
          tc_e7_taMN)
      t_lhs_vielbein = (tf.eye(56, dtype=tf.complex128) +
                        t_d_gen_e7 + 0.5 * t_d_gen_e7 @ t_d_gen_e7)
      t_potential = (
          tf_sugra_tensors(t_v70,
                           tf.constant(True),
                           t_lhs_vielbein,
                           tc_rhs_E))[-1]
    return tape.gradient(t_potential, t_d_v70_onb)
  #
  t_d_v70_onb = tf.Variable(numpy.zeros(70), dtype=tf.float64)
  tape = tf.GradientTape(persistent=True)
  with tape:
    tape.watch(t_d_v70_onb)
    grad_potential = tf.unstack(tf_grad_potential_lhs_onb(t_d_v70_onb))

  t_mm = tf.stack([tape.gradient(grad_potential[k], t_d_v70_onb)
                  for k in range(70)], axis=1)
  stensors = so8_sugra_tensors(t_v70, tc_rhs_E)
  return (t_mm * (36.0 / tf.abs(stensors.potential))).numpy()


### Scanning

def scanner(
    use_dwn_stationarity=True,
    so_pq_p=8,
    seed=1,
    scale=0.15,
    stationarity_threshold=1e-4,
    relu_coordinate_threshold=3.0,
    gtol=1e-4,
    f_squashed=tf.math.asinh):
  """Scans for critical points in the scalar potential.

  Args:
    use_dwn_stationarity: Whether to use the explicit stationarity condition
      from `dwn_stationarity`.
    so_pq_p: SO(p, 8-p) non-compact form of the gauge group to use.
    seed: Random number generator seed for generating starting points.
    scale: Scale for normal-distributed search starting point coordinates.
    stationarity_threshold: Upper bound on permissible post-optimization
      stationarity for a solution to be considered good.
    relu_coordinate_threshold: Threshold for any coordinate-value at which
      a ReLU-term kicks in, in order to move coordinates back to near zero.
      (This is relevant for noncompact gaugings with flat directions,
      where solutions can move 'very far out'.)
    gtol: `gtol` parameter for scipy.optimize.fmin_bfgs.
    f_squashed: Squashing-function for stationarity.
      Should be approximately linear near zero, monotonic, and not growing
      faster than logarithmic.
  Yields:
    `Solution` numerical solutions.
  """
  # Use a seeded random number generator for better reproducibility
  # (but note that scipy's optimizers may themselves use independent
  # and not-easily-controllable random state).
  rng = numpy.random.RandomState(seed=seed)
  def get_x0():
    return rng.normal(scale=scale, size=70)
  #
  tc_rhs_E = tf.constant(get_so_pq_E(so_pq_p), dtype=tf.complex128)
  def f_potential(scalars):
    return so8_sugra_tensors(tf.constant(scalars), tc_rhs_E).potential.numpy()
  #
  f_grad_pot_sq_stationarity = (
      None if use_dwn_stationarity
      else get_tf_stationarity(
          lambda t_pos: so8_sugra_tensors(t_pos, tc_rhs_E).potential))
  #
  def f_t_stationarity(t_pos):
    if use_dwn_stationarity:
      stensors = so8_sugra_tensors(t_pos, tc_rhs_E)
      stationarity = dwn_stationarity(stensors.a1, stensors.a2)
    else:
      stationarity = f_grad_pot_sq_stationarity(t_pos)
    eff_stationarity = stationarity + tf.reduce_sum(
        tf.nn.relu(abs(t_pos) - relu_coordinate_threshold))
    return eff_stationarity
  #
  def f_opt(pos):
    t_pos = tf.constant(pos)
    t_stationarity = f_squashed(f_t_stationarity(t_pos))
    return t_stationarity.numpy()
  #
  def fprime_opt(pos):
    t_pos = tf.constant(pos)
    tape = tf.GradientTape()
    with tape:
      tape.watch(t_pos)
      t_stationarity = f_squashed(f_t_stationarity(t_pos))
    t_grad_opt = tape.gradient(t_stationarity, t_pos)
    return t_grad_opt.numpy()
  #
  while True:
    opt = scipy.optimize.fmin_bfgs(
        f_opt, get_x0(), fprime=fprime_opt, gtol=gtol, maxiter=10**4, disp=0)
    opt_pot = f_potential(opt)
    opt_stat = f_opt(opt)
    if numpy.isnan(opt_pot) or not opt_stat < stationarity_threshold:
      continue  # Optimization ran into a bad solution.
    solution = Solution(potential=opt_pot,
                        stationarity=opt_stat,
                        pos=opt)
    yield solution

### Demo.

def demo(seed=0,
         scale=0.2,
         use_dwn_stationarity=True,
         so_pq_p=8,
         num_solutions=5,
         f_squashed=tf.math.asinh):
  solutions_iter = scanner(scale=scale, seed=seed,
                           use_dwn_stationarity=use_dwn_stationarity,
                           so_pq_p=so_pq_p, f_squashed=f_squashed)
  for num_solution in range(num_solutions):
    sol = next(solutions_iter)
    print('=== Solution ===')
    pprint.pprint(sol)
    mm0 = so8_sugra_scalar_masses(sol.pos, so_pq_p)
    print('\nScalar Masses for: V/g^2=%s:' % sol.potential)
    print(sorted(collections.Counter(
        numpy.round(numpy.linalg.eigh(mm0)[0], 3)).items()))

demo()
